Reconstructing signals from noisy data with unknown signal and noise covariance 
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We derive a method to reconstruct Gaussian signals from linear measurements with Gaussian 
noise. This new algorithm is intended for applications in astrophysics and other sciences. The start- 
ing point of our considerations is the principle of minimum Gibbs free energy which was previously 
used to derive a signal reconstruction algorithm handling uncertainties in the signal covariance. We 
extend this algorithm to simultaneously uncertain noise and signal covariances using the same prin- 
ciples in the derivation. The resulting equations are general enough to be applied in many different 
contexts. We demonstrate the performance of the algorithm by applying it to specific example 
situations and compare it to algorithms not allowing for uncertainties in the noise covariance. The 
results show that the method we suggest performs very well under a variety of circumstances and is 
indeed qualitatively superior to the other methods in cases where uncertainty in the noise covariance 
is present. 
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I. INTRODUCTION 

The problem of signal inference consists of reconstruct- 
ing a set of parameters or even a continuous field s from 
some data set <i, which is influenced in some way by the 
signal, 



d=/(s)+n. 



(1) 



Two problems will arise. First, the function / may not 
be invertible and, second, the noise term n will not be 
known. In the Bayesian framework, one uses prior in- 
formation on the signal and the noise term to calculate 
a best estimate for the true signal realization or, ideally, 
the whole probability distribution for the signal given the 
prior information and the information contained in the 
data. 

Symmetry considerations and knowledge about the un- 
derlying physics of the signal and the measurement pro- 
cess may restrict the class of priors that one has to con- 
sider. They might, however, still contain some free pa- 
rameters that then become part of the inference problem. 
The case in which the signal covariance contains uncer- 
tain parameters was tackled in pQ, producing a whole 
class of filters for this problem. The filter that we extend 
in this work was reproduced in [3], where the principle 
of minimum Gibbs free energy was introduced (cf. also 
Sect. |III| ), and successfully applied in an astrophysical 
setting in |5]. 

Here, we focus on the case where we can assume zero- 
mean Gaussian priors both for the signal and for the 
noise. The priors are therefore completely characterized 
by the respective covariance matrices. Our goal is to ex- 
tend the study of [3J to the case in which both the signal 
covariance and the noise covariance contain parameters 



* niels@mpa-garching.mpg.de 



that are not known a priori. This is motivated mainly by 
applications from the field of astrophysics. The theory 
and resulting filter formulas, however, are of general ap- 
plicability. Gaussian noise, e.g., is omnipresent in nearly 
every area of the natural sciences and the situation in 
which its variance is not precisely known should be a 
rather common one. 

Previous work dealing with the problem of unknown 
noise variance has mainly dealt with specific applications. 
One of these applications is the field of image reconstruc- 
tion. Here, it is usually assumed that the measured pic- 
ture is the sum of the underlying signal and a white Gaus- 
sian noise term. Often, it is further assumed that the 
noise level, i.e. its variance, is the same in every image 
pixel. A comparison of different algorithms for noise esti- 
mation under these assumptions can be found e.g. in [7]. 
An example for an algorithm allowing for inhomogeneous 
noise is presented in [12], where a wavelet transform of 
the image is applied and the lack of correlated noise is 
exploited. Most of these algorithms, however, are not 
derived by rigorous statistical calculations but rather by 
a combination of intuition and experience. 

From a mathematical viewpoint, the problem of an 
unknown noise prior has received some attention in the 
theory of density deconvolution, which deals with the in- 
ference of the probability density for a signal from mea- 
surements with additive noise. Here, the signal is usually 
assumed to consist of independent identically distributed 
variables. The case of Gaussian noise with unknown vari- 
ance has been considered e.g. by [5] and [TT] . 

In this work, we create a general setting with well de- 
fined assumptions and a traceable derivation of a general 
filter formula within a Bayesian framework, not loosing 
sight of its applicability. Our result can accomodate a 
host of different assumptions and models, such as cor- 
related or uncorrelated noise. It allows for a distinction 
between the data space and the signal space, with possi- 
bly different numbers of degrees of freedom. 



The remainder of the paper is organized as follows. In 
Sect.[ll]we introduce our model for the measurement pro- 
cess and the notation that is required. The derivation of 
the filter formulas follows in Sect. IIIII We then demon- 
strate the usefulness of our filter by applying it in a set 
of mock observational situations in Sect. HVl and discuss 
the implications in Sect. [VI 



II. SIGNAL MODEL AND NOTATION 

We assume a linear measurement process where the 
data are a superposition of a linear signal response and 
a noise term, 



d = Rs 



(2) 



Here, both the data and noise and the signal can be finite- 
dimensional vectors or continuous fields defined on some 
manifold. The response matrix R maps a field in the 
signal space to a field in the data space. In the continuous 
limit the matrix vector product becomes 



(R*)i 



Q.X rLj iT Sqr 



(3) 



where the index denotes the value of a field at this posi- 
tion. In physical applications the data vector will always 
be discrete, since only a finite number of measurements 
can be taken, while the signal space might well be contin- 
uous. The result of any numerical signal reconstruction, 
however, will at best be a discretized version of the con- 
tinuous field. 

Further, we assume Gaussian prior statistics both for 
the signal and for the noise contribution, i.e. s ^^ Q(s,S), 
n <-^ Q(n, TV), where 



Q(a,A) 



1 



VWa 



exp 



1 



a*A- 



(4) 



denotes a multivariate Gaussian distribution in a with 
covariance matrix A. We use the dagger symbol to indi- 
cate a scalar product, 



a ] b- 



dx a x b x , 



(5) 



This 



and the asterisk to denote complex conjugation, 
corresponds to the notation introduced in [2]. 

The problem of signal reconstruction is to find an op- 
timal estimate m for the signal realization that the mea- 
sured data arose from. Optimality in an L 2 -norm sense 
leads to 



m 



(s) 



■P(s\d) 



I 



Vs sT(s\d), 



(6) 



i.e. the posterior mean. The integration is performed 
over all possible signal configurations. In the discrete 
case this becomes a product of one-dimensional integrals, 



Vs 



+oo 



dsi 



dso 



(7) 



where s = (si, S2, . • . ) is the vector of signal values at lo- 
cations 1,2, Ideally, we would also like to obtain 

some information on the posterior distribution V(s\d) 
other than its mean. If the signal and noise covariances 
are known, the posterior is a Gaussian Q{s — m, D) with 
mean 

m = Dj, (8) 

and covariance D, where j = R^N~ x d is called the in- 
formation source and D = (S -1 + R^N _1 R) the in- 
formation propagator [cf. 2 and the dagger attached to 
a matrix denotes its hermit ian conjugate. 

In this paper we are concerned with the case in which 
neither the signal covaricance matrix S nor the noise co- 
variance matrix N are known. We parameterize these 
matrices as sums of their eigenvalues pk and fjj multi- 
plied with the projectors onto the respective eigenspaces 
Sk and Nj . The parameters can be rescaled by including 
some numerical values Sk and hj in the project ion- like 
matrices, making the rescaled version of the parameteri- 
zation 



k 



where 



and 



Pk 



Vk_ 

Sk 



Vj 



Sk — s~kSk, Nj 






hjNj. 



(9) 
(10) 

(11) 
(12) 



Furthermore, we define the pseudo-inverse matrices 



n^Nj, 



so that S k 1 Sj £ 



S- 1 = s^Sk and TV" 1 

and N^Nj are identity operators on the respective 

eigenspaces. 

We assume here that the eigenspaces corresponding 
to the different eigenvalues are known a priori, e.g. from 
symmetry considerations. However, the formalism allows 
for eigenvalues of different eigenspaces becoming equal a 
posteriori. 

Finally, we also need to define some priors for the pa- 
rameters pk and r]j . As was done in pQ and [3] , we assume 
each parameter to be a priori independent from all the 
others and use inverse Gamma distributions, i.e. power 
laws with exponential cutoff, as priors for the individual 
parameters, 



i 



v( P ) = n 
v( V ) = n 



q k T(a k - 1) \q k 



Pk 



rjTiPj-l) \r 



V 3 



exp 



exp 



_<lk_ 
Pk 

r A 
Vj 



(13) 
,(14) 

.(15) 



The parameters a k and f3j determine the steepness of 
the power law and the parameters q k and Tj give the po- 
sition of the cutoff. In the limit (a kl j3j) — >• (1,1) and 
(qk^rj) — >> (0,0), this turns into the so-called Jeffreys 
prior, which is flat on a logarithmic scale and can there- 
fore be characterized as non-informative. 



weight beteween the internal energy term and the entropy 
term in the Gibbs free energy. 

Approximating the posterior with a Gaussian with 
mean m and covariance D, 



V(s\d) nQ(s-m,D), 



(20) 



III. DERIVATION OF THE FILTER FORMULAS 

With the priors for s, n, p, and 77, we can calculate the 
joint probability of the signal and the data by marginal- 
izing over the parameters p and 77, 

P(s, d) = Jvpjvr] P(s, d\p, rj)Pfa rj) 

= Jvpjvr] P(d\s^rj)P(8\p)P{p,rj) (16) 

= jvpjvr ] g{d-Rs,N)g{s,S)V{p,r ] ). 

Solving the integrals yields 

r(7fc)C _1 



p(8,d)=n, 



n 



k T(a k -l)(2n) pk/2 
( 1 \ ~ lk 

[Qk + 2 S * S * ls ) 

T(/3 j -l)(27rf j/2 
r 3 + )-(d-Rs) ] N J ~ 1 (d-Rs) 



(17) 



where p k = tr (S k 1 S k ) : Pj = tr (Nj 1 N j ), >y k = p k /2 + 
a k — 1, and £j = /jLj/2-\-/3j — 1. Note that the posterior is 
proportional to this joint likelihood for any given dataset. 
One could construct the maximum a posteriori esti- 
mator, however, this was shown in pQ to perform poorly 
due to a perception threshold, i.e. modes with too little 
power in the data are completely filtered out. A bet- 
ter estimate for the posterior mean of the signal can be 
constructed using the formalism of minimum Gibbs free 
energy, derived in [3 , where thermodynamic quantities 
are introduced by identifying the posterior probability 
density with a canonical density funciton according to 



P(s\d) 



P(s,d) 



P(d) 

The Gibbs energy is then 



-T- 1 H(s y d) 

Z(d) ' 



G = U-TS B , 



(18) 



(19) 



where U = {H)-pr s \ d \ is the internal energy, Sb = 
(—\ogP{s\d))y,, s \ d \ the Boltzmann entropy, and H = 
— logP(s,d) is called the information Hamiltonian. The 
temperature T serves as a tuning parameter, shifting the 



gives an approximate internal energy U, an approximate 
entropy Sb, and therefore an approximate Gibbs energy 

G{m,D) = U(m,D) -TS B 

= (H(s, d)) Q{s _^ D) - |tr (1 + log (2ttD)) . 

(21) 

For T = 1, this approximate energy is, apart from an ad- 
ditive constant, identical to the non-symmetric Kullback- 
Leibler distance [6\ between the full posterior and the 
Gaussian approximation, 



G(m, D) = (H(s, d) + log (G(8 - m, D))) g{s _ mD) 

' g{s — m,D) 



Vs Q(s — ra, D) log 
Vs Q(s — m,D) log 



V(s,d) 
Q(s - m,D) 



V(s\d) 
+ log(V(d)) 

= d KL [G(s - m, D), T(s\d)} + log (Z(d)) , 

(22) 

as was shown already in [3 . 

The approximate internal energy in our case, calcu- 
lated from the joint probability of Eq. (17), is 



*7(m,£)^^7 fc (log 



Qk 



s^S: 



G(s-m,D) 



= -A k 



J2 Sj (log L + \ (d - Rs)l N- 1 (d - Rs) 



Q(s-m,D) 



(23) 



where we have dropped terms that are independent of m 
and D. The logarithms can be expanded in an asymp- 
totic power series, giving 



Ak =\og(q k ) 



(-!)• 



-E^(U * iJs^s qi 



Q(s-m,D) 



■:Aki 



(24) 



and 



B j =log(r j )-J2 



i-iy 



i lr i 

2=1 J 



rj + - (d - Rs) J N~' (d - Rs) - fj 



Q(s-m,D) 



:BA 



(25) 



where we have chosen the linear dependencies to be cap- 
tured by 



and 



qk 



Qk 



|.'V. 



g(s-m,D) 



% + -tr((mm t + Z>)5^ 1 ) 



(26) 



multaneously are 

m = Dj 

qk + |tr ((n 



Pfc 



% 



■jw; 



2 



+ Qfc - 1 



(31a) 
(31b) 



+ |tr (((d - iJra) (d - #m) f + iJDijt) JVT 1 ) 



2 



(31c) 



Thus, we find both the posterior mean and the posterior 
covariance for the signal. Note that the first two of these 
three equations were already found in jj] and [3], where 
the reconstruction of signals with unknown power spec- 
tra is discussed. The term critical filter was coined in pQ 
to refer to this filter since it belongs to a family of filters 
lying on a line in the parameter plane of [1] that separates 
the filters with a perception threshold from those with- 
out. The additional uncertainty in the noise covariance 
that we introduce here simply adds one more equation, 
leading to an extended critical filter. 



-h-(d-iis) t Nr 1 (d-Rs) 

2 / Q(s-m,D) 

~tr (((d - Rm) (d - Rm) ] + RDR^ JVT 1 ) , 

(27) 



respectively. 

Here, we restrict ourselves to the zeroth order solution, 
i.e. we neglect all contributions from A and B. Further- 
more we set T = 1. The case with T ^ 1 is discussed up 
to second order in Appendix |A| 

Now we search for the optimal Gaussian approximation 
to the posterior by minimizing the approximate Gibbs 
energy, which is equivalent to minimizing the Kullback- 
Leibler distance between the two probability densities, 
according to Eq. (22). Taking the functional derivatives 



of Eq. (21) with respect to m and D and equating them 



to zero yields the equations 



m = Dj, 



= £^v^ 



D = Ef^+ES^ 1 * 



qk 



(28) 
(29) 

(30) 



By comparing these expressions to the Wiener filter for- 
mula, Eq. ([8|, we can read off the parameters Pk = — 

and rjj = J 2 - for the signal and noise covariance matrix, 
respectively. 

So altogether the equations that need to be solved si- 



IV. APPLICATION TO SIMULATED SIGNALS 

Here we demonstrate the performance of our signal 
reconstruction algorithm under different circumstances. 

A. Setup 

We consider two different scenarios. First, we consider 
a simple one-dimensional test case, where the signal is 
supposed to be a real field defined over some interval 
with periodic boundary conditions. We discretize the in- 
terval into 2048 pixels. For simplicity, we set the response 
matrix R to be the identity operator, so that our data 
set consists of 2048 individual points as well. We fur- 
ther assume statistical homogeneity for the signal field, 
leading to a covariance matrix that is diagonal in Fourier 
representation, 



Skk' = ( s /e 5 fc')-p( s ) = 3kk'Ps(k), 



(32) 



with the power spectrum P s (k) on its diagonal. For this 
power spectrum we choose a simple power law 



P s (k) oc (1 + fc)" 



(33) 



and draw a random realization of the signal from it. 

Motivated by astrophysical applications, we also con- 
sider a real signal field on the sphere, 



s: S z 



R. 



(34) 



Using again R = 1, the data and noise are also fields on 
the sphere, 



d, n : S 



R. 



(35) 




FIG. 1. Comparison of different filter algorithms in the one-dimensional test case. Each column corresponds to a different 
setting. The signal, drawn from a power law power spectrum, is the same in each case and depicted in each panel with a solid 
line. The left column contains homogeneous noise, while in the middle column, the noise is suppressed in the left third of the 
interval and enhanced in the right third, and in the right column the noise is enhanced in some individual pixels. The first 
row shows the signal realization along with the data. The second row shows the reconstruction using the Wiener filter formula, 
assuming the correct power spectrum and under the assumption of homogeneous noise; the third row shows the critical filter 
reconstruction, assuming the power spectrum to be unknown, but still assuming homogeneous noise. The last row, finally, 
shows the extended critical filter reconstruction in which both the signal power spectrum as well as the noise variance are 
assumed to be unknown. The respective reconstructions are depicted by a dashed line which lies on top of the solid one in 
many cases. In the two right panels of the first row, some of the data points lie outside the area that is shown. 
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FIG. 2. Comparison of different filter algorithms in the spherical case. Each column corresponds to a different setting. The 
signal, drawn from a power law power spectrum, is the same in each case. The left column contains homogeneous noise, while 
in the middle column, the noise is suppressed in the southern third of the sphere and enhanced in the northern third, and in the 
right column the noise is enhanced in some individual pixels. The first row shows the signal realization and the second row the 
data. The third row shows the reconstruction using the Wiener filter formula, assuming the correct power spectrum and under 
the assumption of homogeneous noise; the fourth row shows the critical filter reconstruction, assuming the power spectrum to 
be unknown, but still assuming homogeneous noise. The last row, finally, shows the extended critical filter reconstruction in 
which both the signal power spectrum as well as the noise variance are assumed to be unknown. 



In the numerical implementation, we use the HEALPdQ 
discretization scheme at a resolution of 7V S id e = 16, which 
leads to 3 072 pixels. Under the assumption of statistical 
homogeneity and isotropy, the signal covariance matrix 
in this case becomes diagonal in the basis given by the 
spherical harmonics components, 

S(£m)(£'m') = ( s lmS£>m' )v(s) ~ ^W^mm'C^ (36) 

where Ci are the angular power spectrum components. 
We draw our signal realization again from a power law 
spectrum, 



C £ = (1+£Y 



(37) 



We assume the noise to be uncorrelated in the position 
basis, making the noise covariance matrix diagonal in this 
basis, 



AT- 

1 v nr 



(%%')p( n ) 



5-, a? 

J nn u fii 



(38) 



where h and ft' denote positions on the sphere or on the 
interval, respectively. Within this framework, we con- 
sider three cases for the noise statistics. In the first one, 
we use homogeneous noise with variance cr? = 1/4 in- 
dependent of h. For the second case we divide the data 
space into three zones. In the left/southern third, we 
suppress the noise variance by a factor of nine and in the 
right /northern third we enhance it by a factor of nine, 
while we leave it unchanged in the middle. Finally, in 
the third case we again assume homogeneous noise with 
variance 1/4, but we enhance the variance in five per- 
cent of the pixels, randomly selected, by a factor of 100. 
Both the signal and the three resulting data realizations 
are shown in Fig. [I] for the one-dimensional case and in 
Fig. [2] for the spherical case, along with the results of 
different reconstructions that we discuss next. 



B. Reconstructions 

We first apply the standard Wiener filter formula, the 
results of which are shown in the second row of Fig. [I] for 
the one-dimensional case and in the middle row of Fig. [2] 
for the spherical case. For this we assume the correct 
power spectrum to be known, but we assume homoge- 
neous noise with variance 1/4 in all three cases. In the 
case where this assumption is correct, the reconstruction 
is known to be optimal and this is confirmed by visual 
inspection of the outcome. In the cases with inhomoge- 
neous noise, the Wiener filter fails to completely filter out 
the noise structures in the data in the regions where the 
noise is underestimated and therefore reproduces some of 
them in the reconstruction, as one would expect. This is 



true for the right (northern) third in the middle column 
of Fig. IT] (Fig. |2|, as well as the noisy pixels in the right 
column. The opposite should happen in the left (south- 
ern) third in the middle column of Fig. fl] (Fig. 2l), where 
the noise is overestimated. One would expect that struc- 
tures in the data that are actually due to the signal get 
filtered out. This is actually happening, although it is 
barely visible in the resulting plots. 

Next we assume that the power spectrum is not known 
a priori, i.e. we apply the critical filter. The resulting 
plots are shown in the third row of Fig. [I] for the one- 
dimensional case and the fourth row of Fig. [2] for the 
two-dimensional case. In the one-dimensional case we 
define the Sk operators of Eq. ([9| to be projections onto 
bins of width A^ = 2 in Fourier space, effectively assum- 
ing that the two scales that enter the bin have the same 
power. This binned power is then represented by the 
parameer pk in Eq. ^ . This binning is necessary in one 
dimension since each individual Fourier component con- 
tains only two degrees of freedom jj In the spherical case, 
we can directly use the angular power spectrum compo- 
nents Cg as parameters and the projection-like operators 
Sk become actual projections onto the ^-th angular scale, 
which contains 2^ + 1 degrees of freedom. In both cases 
we assume Jeffreys prior for the unknown parameters. 
Then we simply iterate the first two lines of Eq. (31), 



while keeping the assumption of homogeneous noise with 
variance 1/4. 

In the cases where our assumptions about the noise are 
true, the resulting map is very close to the Wiener fil- 
ter reconstruction, confirming the assessment of [TJ |3j [8] 
that the critical filter can yield a very accurate recon- 
struction, even if the power spectrum is completely un- 
known. In the cases where we have made false assump- 
tions about the noise, however, we see the same prob- 
lems that the Wiener filter reconstruction has, only much 
stronger pronounced. This is because the reconstructed 
power spectrum now actually accounts for the features 
in the data that are due to noise where this is underesti- 
mated. With this power spectrum, the map reconstruc- 
tion tends to favor these features even more than when 
the correct power spectrum is used. This amplifying ef- 
fect is again much more prominent where the noise was 
underestimated than where it was overestimated. 

Finally, we account for the possibility that we might 
have misestimated the noise statistics by applying the full 
extended critical filter, derived in Sect. |III| As projection- 
like matrices Nj we choose projections onto the j-th pixel 
of the interval and sphere, respectively, multiplied with 
our original guess for the noise variance in that pixel, 
a 1 - — 1/4. In this way, the parameters rjj become cor- 
rection factors for the noise variance of each data point. 
For the prior parameters we choose j3j = 2 and we adapt 



1 The HEALPix package is available from http://healpix.jpl. 
nasa.gov 



2 Note that since our signal is real, the Fourier components asso- 
ciated with negative fc-values do not contain additional degrees 
of freedom. 



Vj such that (log {r]j)) v t \ = 0. After iterating the full 
set of equations ( [31] ), we obtain the results shown in the 
bottom rows of Fig. [I] and Fig. [2] In the case with ho- 
mogeneous noise, we still get a result that is similar to 
the Wiener filter one. This shows that we do not lose 
much by allowing for some uncertainty in the noise co- 
variance. In the cases in which our original noise estimate 
was wrong, however, we obtain reconstructions of a much 
higher quality than from the critical filter. Obviously, our 
algorithm succeeds in uncovering the false error bars in 
our dataset and correcting them. This works especially 
well in the case where only individual pixels have under- 
estimated noise variance. This setting makes it especially 
easy for the algorithm to infer the signal statistics from 
all the other pixels and find the pixels in which the data 
points and the signal are inconsistent with one another. 
However, even in the case where one third of the space 
is covered with underestimated noise, our algorithm still 
does a good job in reconstructing the original signal. In 
the spherical case, the extended critical filter performs 
even better than the Wiener filter. This is also true 
for the one-dimensional case in the scenario where the 
noise in individual pixels is enhanced. In the scenario 
with enhanced and suppressed noise in one third of the 
one-dimensional interval, the Wiener filter performs bet- 
ter than the extended critical filter. It should be noted, 
however, that using the Wiener filter is not an option if 
the power spectrum of the signal is not known a priori. 

Some further insight can be gained by looking at the re- 
constructed angular power spectra for the spherical case. 
These are shown in Fig. [3j In the case with homogeneous 
noise, the critical filter recovers the true power spectrum 
almost perfectly, while the extended critical filter misses 
some power on the smallest scales, i.e. some of the small- 
scale power in the data is falsely attributed to noise and 
therefore not represented in the signal power spectrum. 
This effect is small, however, and does not greatly influ- 
ence the resulting map. 

In the case in which the noise is highly inhomogeneous, 
being higher and lower in one third of the data space each, 
the extended critical filter misses quite a lot of power on 
small scales. This results in the slightly oversmoothed 
map seen in Fig. [2] The critical filter, however, that 
operates under wrong assumptions for the noise statis- 
tics, overestimates the power on small scales significantly. 
This is in agreement with the very noisy reconstructed 
map. 

It is in the third case, in which the noise is greatly 
enhanced in individual pixels, that the extended criti- 
cal filter shows its full strength. While the critical filter 
attributes the power in the faulty pixels to the signal 
and therefore overestimates the signal power by orders of 
magnitude, the extended version accounts for the mises- 
timated error bars and does not account for these pixels 
in the signal power spectrum. While the result is a power 
spectrum that is slightly underestimated on the smallest 
scales, this is again only a comparatively small error, still 
leading to a good reconstruction. 
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FIG. 3. Comparison of the different reconstructed angular 
power spectra for the spherical scenario. The solid line de- 
picts the theoretical power spectrum which is also used in the 
Wiener filter reconstructions. The dashed line corresponds to 
the power of the specific signal realization and the dash-dotted 
and dotted lines to the power spectrum reconstructed with 
the critical filter and the extended critical filter, respectively. 
Panel (a) shows the case with homogeneous noise, panel (b) 
the one in which the noise is enhanced and suppressed in one 
third of the sphere each, and panel (c) the one in which the 
noise is enhanced in individual pixels. 



These findings are confirmed by Fig. [4j which shows 
the differences of the nine reconstructions and the sig- 
nal realization. Our extension of the critical filter clearly 
brings the strongest improvement in the case where the 
noise is enhanced in individual pixels, while also lowering 
the error in the case with an extended region of under- 
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FIG. 4. Absolute value of the pixelwise difference between the reconstructed maps and the signal realization for the spherical 
scenario. Each row shows the results for a different filter algorithm. As in Fig. [2] the left column shows the case with 
homogeneous noise, the middle column the one with enhanced noise in the northern third of the sphere and suppressed noise in 
the southern third, and the right column shows the case where the noise is enhanced in individual pixels. Note that the color 
bar differs from the one used in Fig. [2] 



extended critical filter 




0.15 



0.25 



0.35 



0.45 



0.55 



FIG. 5. Pixelwise uncertainty of the extended critical filter reconstructions in the spherical scenario. The left panel shows the 
case with homogeneous noise, the middle panel the case with enhanced noise in the northern third and suppressed noise in the 
southern third, and the right panel the case with enhanced noise in individual pixels. 



estimated noise. The same can directly be seen for the 
one-dimensional case in Fig. [T] 

Finally, we plot the standard deviation per pixel of the 
Gaussian approximation (20) to the posterior probabil- 



ity distribution, i.e. the square root of the diagonal of 



the covariance matrix in the pixel basis, a = ^diag(.D), 
in Fig. [5] for the extended critical filter. This can be 
interpreted as an estimate for the la-error bar of the re- 
constructed maps. The region with enhanced noise in the 
second scenario is clearly marked out by a higher uncer- 
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tainty of the reconstruction due to the corrected entries 
of the noise covariance matrix entering the information 
propagator D. Note, however, that the full posterior is 
non- Gaussian and the la range can therefore not neces- 
sarily be interpreted as a 68% confidence interval, espe- 
cially since we are using a zeroth order approximation 
to calculate the Gaussian approximation. In fact, in our 
spherical example only about 50% of the pixels of the 
reconstructions lie within la of the correct signal in all 
three noise scenarios. 

Fig. [6] extends this study to the one-dimensional case. 
We plot the difference of the signal and the reconstruc- 
tion result of the extended critical filter, along with lines 
depicting ±cr, for the three different noise settings. In 
the case where the noise variance is constant within each 
third of the interval, the one-sigma curve exhibits a step- 
like behavior at the boundaries of the thirds, although 
this effect is relatively small. In the case of homogeneous 
noise, the one-sigma curve is roughly constant while the 
individual noisy pixels in the last scenario are reflected 
in the one-sigma curve by its large variations from pixel 
to pixel. The fraction of the pixels for which the recon- 
struction lies within la of the correct signal in the one- 
dimensional case is 50% in the case with homogeneous 
noise, 63% in the case with enhanced and suppressed 
noise in a third of the interval each, and 69% in the case 
with enhanced noise in individual pixels. 



V. DISCUSSION 

Using the formalism of minimum Gibbs free energy we 
have extended the critical filter algorithm, developed in 
PQ and [3], to an algorithm that allows for uncertainties 
both in the signal covariance and in the noise covariance. 
We have demonstrated the performance of our algorithm, 
Eqs. (31), by applying it to a set of mock observations 



on the sphere, as well as in a simple one-dimensional test 
case. 

These applications have shown that the extended crit- 
ical filter performs outstandingly if only a few individual 
data points have a misestimated error bar. However, even 
in a case where large portions of the data are affected, the 
algorithm was shown to perform inarguably better than 
the critical filter, using a fixed - and faulty - assumption 
about the noise statistics. We have also compared the 
results to those obtained from a Wiener filter reconstruc- 
tion, using the correct power spectrum, which is known 
to be optimal if the assumptions about the noise statis- 
tics are correct. This filter was demonstrated, however, 
to lead to reconstructed maps that are much further from 
the true signal than the results of the extended critical 
filter in some cases where the assumptions are not cor- 
rect. 

The choice of the two-sphere as the space on which our 
signal is defined was motivated by astrophysical applica- 
tions, where we could think of the signal as an all-sky 
field or some quantity defined on the surface of a star or 



a planet. Applications in other fields of physics are abun- 
dant. However, it should be noted that there is nothing 
special about the sphere. We could equally well have cho- 
sen a more-dimensional euclidean space, using the power 
spectrum defined in Fourier space instead of the angular 
power spectrum, as we have done in the one-dimensional 
scenario. 

Furthermore, our choice of the identity operator as re- 
sponse matrix was made only on account of simplicity. 
It allowed us to represent the data in the same fashion 
as the signal. It should be clear, however, that the de- 
rived filter formulas, Eqs. (31 ), are valid for any response 



matrix, even a singular one. Applications of the critical 
filter with non-trivial response matrices were presented in 
[I] and [8] and such a response would not pose a problem 
for the extended version of the filter. 

The problem of signal reconstruction with some un- 
certainty in the noise variance is certainly one of general 
interest. There are several ways in which uncertainty in 
the noise variance might arise. It may be due to ques- 
tionable assumptions that enter in the calculation of the 
error bars of the data. Another possibility is that it arises 
from the definition of the signal itself. The quantity of 
interest may only be part of what has been measured in 
the first place in which case the rest of the data would be 
noise with essentially unknown variance. All these factors 
come together in the reconstruction problem considered 
in [8] . An extension of that reconstruction, using addi- 
tional sets of data and the improved algorithm presented 
here, is planned [9]. 

It should be noted, however, that even with the ex- 
tended critical filter, some knowledge about either the 
parameters in the signal covariance or the ones in the 
noise covariance is needed to arrive at a sensible recon- 
struction. Leaving them both completely free would lead 
to a degeneracy between signal and noise that cannot be 
resolved. Only by assigning an informative prior to at 
least one of the two sets of parameters is this degener- 
acy broken. Furthermore, the functional bases in which 
the signal and noise covariances are diagonal, i.e. their 
eingenspaces, need to differ to allow for a separation of 
noise from signal. 
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FIG. 6. Pixelwise uncertainty of the extended critical filter reconstructions in the one-dimensional scenario. The left panel 
shows the case with homogeneous noise, the middle panel the case with enhanced noise in the right third and suppressed noise 
in the left third, and the right panel the case with enhanced noise in individual pixels. The dark curves represent ±<r and the 
light curve the difference between the reconstruction and the correct signal. 



Appendix A: Higher order solutions 



Here we briefly list the results for nontrivial temper- 
atures up to second order, i.e. considering terms up to 
i = 2 in Eqs. p4|) and (25 ). Since the first order terms are 



zero for our choice of q^ and fj , we list only the resulting 
filter formulas for the zeroth and second order internal 
energy. 



Zeroth order 



The zeroth order solution with arbitrary temperature 
is rather similar to the one with T = 1 presented in 

I 



Sect. Ill It is given by 
m = D'j, 



D' 



Y,f k s ^ + Y.t.^ N i lR 






D = TD'. 



(Al) 
(A2) 

(A3) 

(A4) 



The mean m is completely unchanged. However, the co- 
variance D of the Gaussian approximation is now T times 
the information propagator D' , i.e. the Gaussian approx- 
imation becomes wider at higher temperature. 



2. Second order 



The second order solution is given by 



m = D'j, 



Si 



j = J2 U fY jR ^d, 



(A5) 
(A6) 



D '= El^+E^v 1 * 



Qk 



X k = l + Itr ((mm* + l -DJ S^DS^ - ±-S?D, 
Yj = l+ 4 tr ( ( (d - Rm) (d - Rm)" 1 + ]-RDR) J N^RDR^N'A - -^R^N^RD, 

D = T\ D'- 1 - V ^S^ 1 (mm)) S^ 1 - V %R ] N~ l ((d - Rm) (d - Rm)^) N^R 



(A7) 

(A8) 
(A9) 

(A10) 
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Again, the only effect of the temperature is to broaden 
the approximate Gaussian. However, in the second order 
solution the operators X^ and Yj appear, destroying the 
one-to-one correspondence between the terms in these 
expressions for D, D' , and j and the Wiener filter for- 



mula Eq. ([8|. Therefore, the values of the parameters pk 
and r]j are not immediately determined by these equa- 
tions. Note, however, that the goal was not to determine 
the signal and noise covariance matrices, but to find the 
optimal Gaussian approximation to the signal posterior, 
given by m and D. 
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